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A  COMPUTER  PROGRAM  FOR  THE  CALCULATION 
OF  DAUBECHIES  WAVELETS 

1.  Introduction 

Recently,  it  has  been  found  that  there  exist  other  sets  of  functions  that  can  serve  as  an 
orthonormal  basis  for  the  set  of  all  square  summable  functions,  L~{R)  [1-6].  In  particular, 
it  has  been  shown  that  for  all  g{t)  6  L'^{R),  and  for  a  certain  prescribed  function 

9it)='  E  gj,ky/^rI’{2nt-2-^k)),  (1) 

jMZxZ 

where 

gj,k  =  r  g{tW{2^{t  -  2-^k))dt.  (2) 

J—oo 

and  Z  denotes  the  integers.  Here,  the  basis  set  is  —  2“-'A-))}^,i.g2x2  and  the 

function  t^(t)  is  called  the  wavelet.  Furthermore,  the  function  v{t)  is  significant  only  over  a 
small  (compact)  portion  of  the  real  line.  Therefore,  it  is  recognized  that  this  representation 
lias  a  sense  of  ‘time  locality.’ 

For  i}{t)  to  be  admissible  as  a  basis  function,  its  Fourier  transform  must  obey  certain 
properties.  In  particular,  for  an  admissible  t/’CO  it  can  be  shown  that 

=  Hf)  =  *'({)®(5)'  (3) 

where 

hV)  =  e-‘^->H-(f  +  l/2),  (4) 

mf)=  E  (5) 

Jfc=— oo 

e(/)  =  E //(2->’/),  {c) 

p=i 

and  {h{k)}k^z  is  a  sequence  such  that  the  following  properties  hold; 

(i)  \mo)\  =  h 

(ii)  h{k)  ~  0{k‘^)  as  k  oo, 

(m)|H(/)P  +  |ff(/+l/2)P=  1. 
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(iv)  \H(S)\  5^  0  for  /  e  (0, 1/2). 


Xote  that  the  inverse  Fourier  transform  of  ©(/)  is  called  the  scaling  function,  i.e., 

e{t)  =  r  Q{f)e^^^^df.  (7) 

J  — OO 

The  difference  between  the  wavelet  and  its  associated  scaling  function  is  rooted  in  the 
difference  between  spanning  different  subspaces  that  compose  the  space  L^{R)-  It  can  be 
shown  that  one  can  define  a  series  of  subspaces  Vj  for  j  €  Z  such  that  \Jjez  ^  j  is  dense  in 
L-{R),  Vj  C  Vj+i,  g{t)  6  Vj  if  and  only  if  g{2t)  €  Vj+i  for  all  j  €  Z,  and  g(t)  €  Vj  if  and 
only  if  g(t  —  2~^k)  G  Vj  for  all  j,k  6  Z.  The  significance  of  the  scaling  function  is  that 
the  set  {y/^9{2^{t  —  k))}k^z  spans  the  subspace  Vj.  On  the  other  hand,  one  can  show  that 
there  exists  a  subspace  Oj  composed  of  functions  that  are  orthogonal  to  those  composing 
1  j  such  that 

0,0''',  =  V>1,  (8) 

where  0  denotes  the  Cartesian  product.  Thus,  one  can  show  that 

U  O;  =  LHR)-  (9) 

;€Z 

The  significance  of  the  wavelet  is  that  the  set  {'/^il!{2^{t  —  2~^k))]i;q^z  spans  the  subspace 
Oj.  .More  generally,  the  set  (\/^U'(2^(t  -  spans  L'^{R). 

The  results  presented  above  suggests  that  one  can,  to  some  degree,  control  the  shape 
of  the  wavelet  in  the  time  domain  according  to  how  one  chooses  the  sequence  {li{k)}k^z- 
One  desirable  property  is  to  have  a  wavelet  with  compact  support  in  the  time  domain,  i.e., 
it  is  time  limited  in  that  it  is  nonzero  only  over  a  given  inteiATil.  Such  a  wavelet  gives  a  true 
sense  of  time  locality.  A  set  of  orthogonal  wavelets  with  compaict  support  was  discovered  by 
Daubechies  [2].  They  are  parameterized  by  an  integer  n,  are  real  valued,  and  are  denoted 
.'IS  for  rz  >  2.  In  feict, 

suppi'n  C  [{I  -  n),n].  (10) 

These  wavelets  are  derived  by  choosing  the  sequence  {h{k)}k€Z  so  that  it  is  of  finite  length. 
The  result  is  the  set  of  sequences  for  n  =  1,2,3,  —  Details  of  the  procedure 

for  finding  these  sequences  can  be  found  in  Daubechies’  original  paper  [2]. 

Daubechies  wavelets  posses  other  desirable  properties.  It  can  be  show  that  thcj’  are 
bounded,  continuous  functions  for  all  n,  and  for  n  >  4,  they  are  continuously  differentiable. 
Furthermore,  for  n  >  4,  Daubechies  wavelets  have  a  finite  spectral  spectral  variance,  i.e., 
for  'Pn(/)  then  ^ 

r  /"l't’n(/)|V/  <  OO.  (11) 

J— 130 

Proof  of  these  properties  can  be  found  in  [4]. 
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Basic  Method  of  Calculating  Wavelets 


Equations  (3)  to  (6)  suggest  how  one  can  calculate  the  Fourier  transform  of  a  Daubechies 
wavelet  and  its  associated  scaling  function.  For  ©„(/)  ^n(0>  one  first  calculates  (approx¬ 

imates)  the  Fourier  transform  of  the  scaling  function  via  the  equation 

Qn(/)  «  n  ^n{2-V),  (12) 

p=0 

where  ^  ^ 

Hn(/)  =  E  (13) 

fc=0 

Once  ©„(/)  is  found,  one  calculates  the  Fourier  transform  of  the  orthogonal  wavelet  as 

'I'n(/)'=  A'n(|)e„(|),  (14) 

where 

A'„(/)  =  e-^2'^//;:(/+l/2).  (15) 

The  truncated  product  in  Eq.  (12)  gives  good  results  for  P  =  20  for  low  values  of  n  {n  =  3), 
to  P  =  25  for  high  values  of  n  (n  =  13).  This  was  checked  by  calculating  the  normalized 
cross  correlation  between  two  Daubechies  wavelets  of  order  n,  where  one  was  derived  by 
using  P  =  iV,  and  the  other  with  P  =  iV  -f  1.  For  P  =  25  (or  P  =  20  for  low  values  of  n) 
the  correlation  was  negligibly  different  from  1. 

Once  the  Fourier  transforms  of  the  Daubechies  wavelet  and  scaling  function  have  been 
calculated,  one  can  find  the  associated  time  domain  functions  by  invoking  the  inverse  Fourier 
transform.  This  can  be  accomplished  efficiently  through  the  use  of  a  fast  Fourier  transform 
(FFT). 

The  program  listed  in  Appendix  A  uses  the  approach  outlined  above  to  calculate 
Daubechies  wavelets  and  their  scaling  functions,  and  is  written  in  VAX  extended  FOR¬ 
TRAN.  The  program  produces  four  sequential  ASCII  files  containing  sampled  versions  of 
the  functions  and  fl„(t)  and  the  magnitudes  of  their  Fourier  transforms.  Certainly  one 
can  modify  the  program  to  produce  the  complex  Fourier  transform.  An  ASCII  file  is  also 
produced  containing  the  parameters  input  to  the  program  by  the  user. 


3.  Running  the  Program 

Program  WAVE  is  designed  to  be  run  interactively,  and  produces  sampled  versions 

of  the  Daubechies  wavelets  and  their  associated  scaling  functions  for  7i  =  2,3 . 15.  It  is 

written  using  double  precision,  and  can  take  several  minutes  of  wall  clock  time  to  run.  The 
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program  makes  use  of  a  radix  2  FFT  as  an  efficient  way  of  numerically  calculating  a  Fourier 
integral. 

The  program  needs  only  one  input  file,  COEFF.DAT,  which  is  listed  in  Appendix  B. 
This  file  contains  the  finite  length  sequences  {fin(A;)}t=o^-  All  other  required  parameters  are 
input  from  the  terminal.  They  are: 

1.  n:  The  order  of  the  wavelet,  an  integer  form  2  to  15. 

2.  npower:  The  power  of  2  yielding  the  FFT  size.  For  example,  npower  =  5  implies  that 
the  FFTs  used  in  the  program  are  of  size  2^. 

3.  Time  domain  sampling  rate  in  Hertz  of  and  ^n(0- 

4.  iter\  The  number  of  product  iterations  used  to  calculate  the  Fourier  transform  of  the 
scaling  function  as  given  in  Eq.  (12),  and  is  the  equal  to  P  in  that  equation.  Usually 
a  value  of  20  to  25  is  a  good  choice. 


These  inputs  are  written  to  the  file  W.WE.ECHO.  Therefore,  the  file  serves  as  a  record  of 
a  single  program  run. 

.A.mong  the  remaining  four  ASCII  files  produced  by  the  program  are  WAVELET. TIME 
and  SCALE.TIME,  which  contain  the  sampled  versions  of  the  wavelet  and  scaling  function. 
Specifically,  they  contain  Is)  atid  9^(1/ /,)  where  i  is  an  integer  and  <  i  < 

2npo..r-\  _  1  Similarly,  WAVELET.SPEC  and  SCALE.SPEC  contain  and  Qn{if,) 

respectively,  where  0  <  i  <  2"'”"'"'  —  1.  All  files  contain  ordered  samples  of  the  functions 
in  column  form.  For  example,  WAVELET.TIME  contains  two  columns  of  numbers  in  the 
following  form: 

2  v((l  -  2'’'”“'"-')//,) 

3  tl;((2-2"'^"-’)//,) 

4  t(’((3-2"P^"-’)//,) 


-  1)//,) 


Figures  1  through  4  show  plots  of  the  contents  of  the  four  output  files  generated  by 
WAVE  for  n  =  5,  npower  =  9,  /,  =  30Hz.,  and  iter  =  25. 
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Figure  2.  The  contents  of  the  file  WAVELET.SPEC  which  is  a  sampled  version  of  the 
wavelet  spectrum  Only  the  first  100  points  are  shown. 
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0.0  2S.0  50.0  75.0  100.0 

Figure  4.  The  contents  of  the  file  SCALE. SPEC  which  is  a  sampled  version  of  the  scaling 
function  spectrum  Only  the  first  100  points  are  shown. 
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Appendix  A 


FORTRAN  SOURCE  CODE  FOR  PROGRAM  WAVE 


program  wave 

c  This  program  calculates  a  wavelet  and  its  spectrum.  Note  that 
c  the  sign  (of  the  exponential  argument)  of  the  Fourier  transform 
c  conforms  to  the  conventional  definition  found  in  the  engineering 
c  literature. 

implicit  none 

real*8  h(50) ,whigh.wlow,deltaw,deltaw2, delta, w,fhigh,ts 
real*8  zreal(8192) ,yreal(8192) 
complex*16  htrans.gtrans,x,y(8192) ,z(8192) 
integer  n.npower , iter , nsize , i , j 

external  htrans 
external  gtrans 

c  Open  output  files  containing  the  spectrum  and  time  series  of  the 
c  wavelet  and  scaling  functions. 

open(unit=20,f ile=’scale.spec’ ,status=’new’ , 

4  access=’sequential’) 

open (unit=10 ,f ile= ’wavelet . spec ’ ,status= ’new’ , 

4  access= ’ sequential ’ ) 

open(unit=31,file=’scale.time’ ,status=’new’ , 

4  access=’sequential’) 

open(unit=41 ,f ile=’wavelet .time ’ , status= ’new ’ , 

4  access= ’ sequential ’ ) 

c  Retrieve  impulse  response. 

write(6 , 1000) 

1000  formate/,’  enter  wavelet  order  n:  ’,$) 


read(5, iOlO)  n 
1010  format (bn, iS) 

call  read_coeff (h,n) 

write(6,4000) 

4000  formate/,’  Enter  sampling  frequency:  ’,$) 
read(5,4010)  fhigh 
4010  format(fl0.7) 
ts  =  1.0/fhigh 
whigh  =  2.0*3.141593*fhigh 

write(6,4020) 

4020  formate/,’  Enter  power  of  2  for  number  of  samples;  ’,$) 
read(5,4030)  npower 
4030  f ormat(bn, i5) 

write(6,4040) 

4040  formate/,’  Enter  number  of  iterations:  ’,$) 
read(5,4020)  iter 

if  (npower  .It.  1  .or.  npower  .gt.  13)  then 
type  *,  ’  power  of  2  out  of  range  ’ 
stop 
endif 

c  Open  output  file  and  echo  input  data. 

open(unit=50 ,f ile=’wave . echo ’ ,status= ’new’ , 
b  access= ’ sequential ’ ) 

write (50, 1900)  n, (2**npower) .npower, fhigh, whigh , ts , iter 
1900  formate’  From  progreim  wave:  ’, 

&//,’  wavelet  order  =  ’.i3, 

4  /,’  FFT  size  =  ’,i5,’  (  power  of  2  =  ’,i3.’  )’, 

4  /,’  Maximum  frequency  of  calculated  spectrum  =  ’,ell.5,’  Hz’, 

4  /,’  =  ’,ell.5.’  Rad/S’. 

4  /,’  sampling  rate  =  ’,ell.5,’  S’, 

4  /,’  number  of  iterations  for  calculating  spectrum  =  ’,i4) 
closeeunit=50) 
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c  Calculate  Fourier  transform. 


nsize  =2**npower 
wlow  =0.0 

deltaw  =  (whigh-wlow)/f loatCnsize-l) 
deltaw2  =  delta/2.0 
w  =  wlow 

j  =  1 

write(6,5610) 

5610  formatC/.Sx, '  «  Calculating  scaling  function  spectrum  >>') 
do  i  =  0, (nsize-1) 

w  =  deltaw*float(i)+wlow 
call  phi_spec(x,w,h,n,iter) 
y(i+l)  =  x*(whigh-wlow) 
write(20,5000)  i,cdabs(x) 

5000  f ormat (1 (2x , i6, 2x , ell . 5)) 

enddo 

write(6,5611) 

5611  format(/,5x, ’  <<  Calculating  wavelet  spectrum  >>') 
do  i  =  0, (nsize-1) 

w  =  deltaw*f loat(i)+wlow 
call  phi_spec(x,w/2.0,h,n,iter) 
z(i+l)  =  -1 .0*x«‘gtrans(w/2.0,h,n) 
write(10,5000)  i,cdabs(z(i)) 
enddo 

c  Note  that  the  use  of  the  IFFT  here  is  only  as  a  numerical 
c  integration.  Furthermore,  note  that  you  are  integrating 
c  only  over  the  right  half  of  the  spectrum,  consequently  the  resulting 
c  impulse  response  is  complex  (the  analytic  signal) .  Note  that  the 
c  true  impulse  response  is  equal  (proportional)  to  the  real  part, 
c  The  wavelet  and  scaling  functions  are  normalized  (unit  energy) 
c  prior  to  writing  them  to  the  output  files. 

do  i  =  l,nsize 

y(i)  =  ((-1.0)**(i-l))*y(i) 
z(i)  =  ((-1.0)**(i-l))*z(i) 
enddo 

call  ifft(y .npower) 
call  ifft(z, npower) 
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do  i  =  l.nsize 

yreal(i)  =  dreaKyCi)) 
zreaKi)  =  dreaKzCi)) 
enddo 

call  normalize(yreal,ts ,nsize) 
call  normalizeCzreal ,ts .nsize) 

do  i  =  l.nsize 

writeOl  ,5000)  (i-l-nsize/2)  ,yreal(i) 
write(41 ,5000)  (i-l-nsize/2) ,zreal(i) 
enddo 

end 

subroutine  read.coeff (h,n) 

c  This  subroutine  reads  the  impulse  response  from  the 
c  data  file  coeff.dat. 

integer  i,n 

realms  h(l) 

open(file=’coeff .dat’ ,status=’old’ ,access=’sequential’ , 
*/,  unit=99) 

do  while  (i  .ne.  n) 
read(99,1000)  h(l),i 
1000  formatCf 14.7,i5) 
enddo 

read(99, 1010)  (h(i),  i  =  2,2*n) 

1010  formatCf 14.7) 

write(6,1999) 

1999  formate/,’  spectrum  coefficients:’) 
write(6,2000)  ((i-l),h(i),  i  =  l,2*n) 

2000  format(2x,’h(’,i2,’)  =  ’.flO.7) 

close(unit=99) 

return 

end 
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complex*16  function  htrans(w,h,n) 


integer  n 
real*8  h(l),w 

htrans  =  (0.0, 0.0) 
do  i  =  l,(2>'‘n) 

htrans  =  htrans+h(i) 

*/,  ♦dcmplx(cos(f  loat(i-l)*w)  ,-sin(f  loat  (i-l)*w) ) 

enddo 

htrans  =  htrans/(l . 414214 .0 . 0) 

return 

end 

complex’»16  function  gtrans(w,h,n) 

integer  n 
real*8  h(l) ,w 
complex*16  htrans 

external  htrans 

gtrans  =  dcmplx(cos(w) ,-sin(w))*conjg(htrans((w+3. 141593) ,h,n)) 

return 

end 

subroutine  phi_spec(x, w,h,n , iter) 

c  This  subroutine  calculates  the  spectrum  of  the 
c  wavelet  (as  in  Daubechies'  definition)  at  the 
c  frequency  w. 

integer  n,iter,i 
complex*16  x, htrans 
real»8  h(l) , w 
external  htrans 

X  =  (1.0, 0.0) 
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do  i  =  l.iter 

X  =  x*htrans((w/(2.0**i)) ,h,n) 
enddo 

return 

end 

subroutine  fft(x,in) 

c  This  subroutine  calculates  an  FFT  of  size  2**m.  It  is  an 
c  ‘in-place’  algorithm. 

complex*16  x(l).u,w,t 
n  =  2**m 

pi  3.14159265358979 
do  20  1  =  l.m 
le  =  2**(m+l-l) 
lei  =  le/2 
u  =  (1.0, 0.0) 

w  =  dcmplx(cos(pi/float(lel)) ,-sin(pi/float(lel) ) ) 
do  20  j  =  l.lel 
do  10  i  =  j ,n,le 
ip  =  i  +  lei 
t  =  x(i)  +  x(ip) 
x(ip)  =  (x(i)  -  x(ip))*u 
10  x(i)  =  t 

20  u  =  u*w 

nv2  =  n/2 
nml  =  n-1 
J  =  1 

do  30  i  =  l,nml 

if  (i  .ge.  j)  go  to  25 

t  =  x(j) 

x(j)  =  x(i) 

x(i)  =  t 

25  k  =  nv2 

26  if  (k  .ge.  j)  go  to  30 
j  =  j-k 

k  =  k/2 
go  to  26 
30  j  =  j+k 

return 
end 
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subroutine  ifft(x,n) 


c  This  subroutine  calculates  the  inverse  FFT  of  the  array  x. 

complex*16  x(8192) 
integer  n.nsize.i 
real*8  realn 


nsize  =  2**n 
realn  =  float(nsize) 

do  i  -  1, nsize 

x(i)  -  conjg(x(i)) 
enddo 

call  fft(x,n) 

do  i  =  1, nsize 

x(i)  =  conjg(x(i)/realn) 
enddo 
return 
end 

subroutine  normalizeCx.ts, nsize) 

c  This  subroutine  energy  normalizes  a  real  time  series. 

integer  nsize 
real*8  x(nsize) ,ts ,sum 

sum  =0.0 
do  i  =  1, nsize 

sum  =  sum+x(i)**2 
enddo 

sum  =  sqrt(sum) 
do  i  =  1, nsize 
x(i)  =  x(i)/sum 
enddo 

return 

end 
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Appendix  B 


SEQUENCES  FOR  DAUBECHIES  WAVELETS 


Listed  below  is  the  file  COEFF.DAT  used  by  program  WAVE.  It  contains  the 
sequences  {h„{k)}l'^^  for  n  =  2,3 15.  It  is  formatted  according  to  the  state¬ 
ment  f ormatCf  14 . 7, 15).  The  first  column  lists  the  values  of  h„{k).  The  second 
column  lists  n,  and  marks  the  begining  of  the  sequence.  For  example,  for  n  =  3, 
we  have  h^il)  =  0.332670552950,  /i3(2)  =  0.806891509311,  /;3(3)  =  0.459877502118, 
fi3(4)  =  -0.135011020010,  /i3(5)  =  -0.085441273882,  and  h^iG)  =  0.035226291882. 


.482962913145  2 

.836516303738 
.224143868042 
-.129409522551 
.332670552950  3 

.806891509311 
.459877502118 
-.135011020010 
-.085441273882 
.035226291882 
.230377813309  4 

.714846570553 
.630880767930 
-.027983769417 
-.187034811719 
.030841381836 
.032883011667 
-.010597401785 
.160102397974  5 

.603829269797 
. 724308528438 
.138428145901 
-.242294887066 
-.032244869585 
.077571493840 
-.006241490213 
-.012580751999 
.003335725285 
.111540743350  6 

.494623890398 
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.751133908021 

.315250351709 

.226264693965 

.129766867567 

.097501605587 

.027522865530 

.031582039318 

.000553842201 

.004777257511 

.001077301085 

.077852054085 

.396539319482 

.729132090846 

.469782287405 

. 143906003929 

.224036184994 

.071309219267 

.080612609151 

.038029936935 

.016574541631 

.012550998556 

.000429577973 

.001801640704 

.000353713800 

.054415842243 

.312871590914 

.675630736297 

.585354683654 

.015829105256 

.284015542962 

.000472484574 

. 128747426620 

.017369301002 

.044088253931 

.013981027917 

.008746094047 

.004870352993 

.000391740373 

.000675449406 

.000117476784 

.038077947364 

.243834674613 

.604823123690 


.657288078051 
. 133197385825 
-.293273783279 
-.096840783223 
. 148540749338 
.030725681479 
-.067632829061 
.000250947115 
.022361662124 
-.004723204758 
-.004281503682 
.001847646883 
.000230385764 
-.000251963189 
.000039347320 
.026670057901  10 

. 188176800078 
.527201188932 
.688459039454 
.281172343661 
-.249846424327 
-.195946274377 
,127369340336 
.093057364604 
-.071394147166 
-.029457536822 
.033212674059 
.003606553567 
-.010733175483 
.001395351747 
.001992405295 
-.000685856695 
-.000116466855 
.000093588670 
-.000013264203 
0.018692339500  11 

0.144048360129 
0.449822419238 
0.685506451221 
0.411710892303 
-.162485521339 
-.274320974144 
0.066025638763 
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0.149791844607 

-.046504355457 

-.066445800596 

0.031336714900 

0.020839548328 

-.015365977170 

-.003339972936 

0.004928945867 

-.000308709907 

-.000893056839 

0.000249184997 

0.000054438816 

-.000034637754 

0.000004494745 

0.013114280902  12 

0.109587064387 

0.377449392844 

0.657445006413 

0.516294170295 

-.044313624533 

-.315809615475 

-.023471399498 

0.182806918672 

0.005686977952 

-.096186633657 

0.010995853244 

0.041627451082 

-.012180151045 

-.012829445168 

0  006713258423 

0.002249393038 

-.002179176553 

0.000006459278 

0.000388621871 

-.000088486615 

-.000024241195 

0.000012775434 

-.000001528836 

0.009204916897  13 

0.082889405900 

0.312115898739 

0.611313131287 

0.589096065406 
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0.086639694877 

-.316237370186 

-.126430468961 

0.177816118862 

0.071915527849 

-.106342427892 

-.026758244166 

0.056034390582 

0.002363616024 

-.023833745174 

0.003917927648 

0.007254616037 

-.002760408506 

-.001315670455 

0.000932006061 

0.000049301053 

-.000165090932 

0.000030664729 

0.000010440501 

-.000004699171 

0.000000521846 

0.006547491642 

0.063360170581 

0.259953209778 

0.569486757657 

0.659765991407 

0.253248224211 

-.245883485949 

-.207221475070 

0.141972692112 

0.144030955893 

-.083519992219 

-.071278880702 

0.054864716315 

0.027555092282 

-.029754599557 

-.005754062318 

0.012711190182 

-.000664409841 

-.003831834380 

0.001038385046 

0.000708200880 

-.000381870689 


-.000042656957 

0.000068164760 

-.000010124883 

-.000004370468 

0.000001706613 

-.000000176357 

0.004129483519 

0.041776513906 

0.179077334434 

0.407237141256 

0.472918366318 

0.110408203455 

-.385470269520 

-.386018759764 

0.027435955232 

0.154106559944 

-.072306793965 

-.118765098585 

0.038151205753 

0.052121808971 

-.030285642429 

-.020301136317 

0.016290846193 

0.004432020222 

-.006931177432 

0.000000520974 

0,002017749394 

-.000450763331 

-.000363203984 

0.000171774677 

0.000023620852 

-.000029991608 

0.000003944981 

0.000001876185 

-.000000681954 

0.000000067409 


